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ABSTRACT 

A  three  dimensional  staggered  grid  model,  based  on  the  Boussinesq 
equations,  is  developed,  and  applied  to  mesoscale  waves.  The  model  is 
tested  with  a  three  dimensional  stable  gravity  wave  and  forecasts  for 
40  time  steps  are  compared  with  the  analytical  solutions  and  also  with 
a  full  grid  model.  Comparisons  show  very   good  agreement  with  both 
forecasts  for  the  stated  time.  Failure  of  the  relaxation  scheme  to 
meet  the  convergence  criterion  after  47  time  steps  is  discussed. 


TABLE  OF  CONTENTS 

I.  INTRODUCTION 8 

II.  MODEL  DEVELOPMENT  EQUATIONS  10 

III.  FINITE  DIFFERENCING  13 

IV.  BOUNDARY  CONDITIONS — 15 

V.  INITIAL  CONDITIONS 20 

VI.  RESULTS 22 

VII.  CONCLUSIONS  AND  SUMMARY —  36 

APPENDIX  A  -  EXTRAPOLATED  LIEBMANN  RELAXATION  SCHEME  FOR 

THREE  DIMENSIONAL  APPLICATION  37 

APPENDIX  B  -  DEVELOPMENT  OF  THE  FORCING  FUNCTION  43 

APPENDIX  C  -  COMPLETE  DEVELOPMENT  OF  THE  BOUNDARY  CONDITIONS  45 

LIST  OF  REFERENCES - 48 

INITIAL  DISTRIBUTION  LIST  4g 

FORM  DD  1473 50 


LIST  OF  FIGURES 


1  Illustration  of  the  grid  system  in  the  horizontal 
direction 16 

2  Illustration  of  the  grid  system  in  the  vertical 

direction 18 

3a     U  component  of  stable  gravity  wave  in  x  direction 

for  J  =  5,  K  =  7  — - 24 

3b     V  Component  of  stable  gravity  wave  in  x  direction 

for  J  =  5,  K  =  7 25 

3c     W  component  of  stable  gravity  wave  in  x  direction 

for  J  =  5,  K  =  7 - 26 

3d     9  wave  in  x  direction  for  J  =  5,  K  =  7 27 

4a     U  component  of  stable  gravity  wave  in  y  direction 

for  I  =  5,  K  =  7 28 

4b     V  component  of  stable  gravity  wave  in  y  direction 

for  I  =  5,  K  =  7 —  29 

4c     W  component  of  stable  gravity  wave  in  y  direction 

for  I  =  5,  K  =  7  — 30 

4d     e  wave  in  y  direction  for  I  =  5,  K  =  7 31 

5a     U  component  of  stable  gravity  wave  in  z  direction 

for  I  =  5,  J  =  5 32 

5b     V  component  of  stable  gravity  wave  in  z  direction 

for  I  =  5,  J  =  5 33 

5c     W  component  of  stable  gravity  wave  in  z  direction 

for  I  =  5,  J  =  5 34 

5d     e  wave  in  z  direction  for  I  =  5,  J  =  5 35 

Al     Three  dimensional  schematic  of  the  location  of 

various  fields  38 

p 
A2     y  Operator  in  three  dimensional  form  assuming 

Ax  =  Ay  =  az 42 


LIST  OF  SYMBOLS  AND  ABBREVIATIONS 

c  Phase  speed  of  disturbance  (m  sec"  ). 

c  Specific  heat  of  air  at  constant  pressure  (joule  gm~  °K~  ). 

f  Cori oil's  parameter  (sec"  ). 

-4    -1 

f  Constant  Coriolis  parameter  (.8365  x  10   sec  ). 

F  Forcing  function  for  the  relaxation  equation. 

g  Acceleration  due  to  gravity  (m  sec  ). 

i  Grid  index  in  x  direction 

IM  Number  of  grid  points  in  the  x  direction. 

j  Grid  index  in  y  direction. 

JM  Number  of  grid  points  in  the  y  direction. 

k  Grid  index  in  z  direction. 

KM  Number  of  grid  points  in  the  z  direction. 

X  Wave  number  in  the  x  direction  =  2  tt/L  . 

L  Wave  length  in  the  x  direction. 

L  Wave  length  in  the  y  direction. 

L  Wave  length  in  the  z  direction. 

m  Wave  number  in  the  z  direction  =  tt/L  . 

P  Pressure  (mb). 

PQ  1000  (mb). 

R  Gas  constant  for  air  (joule  gm~  °K~  ). 

S  Used  to  denote  any  arbitrary  scalar  quantity 

t  Time  (sec). 

T  Temperature  (°K). 

u  x  component  of  horizontal  velocity  (m  sec"  ). 

v  y  component  of  horizontal  velocity  (m  sec"  ). 


\T3  Three  dimensional  wind  vector. 

w  Vertical   component  of  velocity  (m  sec"   ). 

W  Scale  value  of  vertical  motion. 

Ax  Mesh  length  in  the  x  direction. 

Ay  Mesh  length  in  the  y  direction. 

Az  Mesh  length  in  the  z  direction. 

C  Vertical   component  of  relative  vorticity  (sec"   ). 

e  Deviation  from  adiabatic  potential   temperature  (°K) 

e0  300  (°K). 

e"  Stratification  potential   temperature   (°K). 

e  Potential   temperature   (°K). 
R/cp. 

u  Wave  number  in  the  y  direction  =  tt/L . 

2       -1 
v  Kinematic  viscosity  coefficient  (m  sec     ). 

$  Pressure  parameter. 

3  Hydrostatic  part  of  the  pressure  parameter. 

3'  Non-hydrostatic  part  of  the  pressure  parameter. 

p  Density  of  air  (kgm  m     ). 

v  Del  operator  =  T|— +j|—  +  F-|-    . 

oX      ay  dZ 

2  3      3      3 

v~     Laplacian  operator  =  — =■  +  — ~  +  — j   • 

3x    3y    3z 

to     Finite  difference  form  of  the  del  operator. 


ACKNOWLEDGMENT 

The  author  wishes  to  express  his  gratitude  to  Dr.  R.  L.  Alberty, 
whose  suggestion  it  was  to  undertake  this  project,  and  whose  unfailing 
optimism,  infinite  patience  and  thorough  understanding  of  the  subject 
of  atmospheric  convection  contributed  greatly  to  the  completion  of  this 
project.  Thanks  also  go  to  Drs.  R.T.  Williams  and  R.L.  Elsberry  for 
discussions  of  dynamics  and  finite  differencing  schemes  used  in  the 
program  and  also  to  Dr.  Williams  for  his  careful  review  of  the  original 
manuscript.  Computer  time  for  the  numerical  calculations  was  provided 
by  the  W.R.  Church  Computer  Center  at  the  Naval  Postgraduate  School. 


I.   INTRODUCTION 

In  the  past  few  years  many  attempts  have  been  made  to  understand 
the  complex  dynamic  interactions  between  the  thunderstorm  and  its 
environment.  Newton  and  Newton  (1959)  proposed  a  dynamical  thunder- 
storm model  which  proposed  that  strong  vertical  motions  within  the 
thunderstorm  cell  create  an  effective  barrier  to  the  flow  of  the 
environmental  air. 

Ogura  (1963)  first  discussed  the  probability  of  using  the  primitive 
equations  for  forecasting  small  scale  dynamical  features.  Interaction 
between  the  convective  elements  and  the  environmental  wind  field  can 
be  investigated  numerically  by  the  use  of  these  equations.  Deardorff 
(1970)  made  a  numerical  study  of  turbulent  flow  using  a  three-dimensional 
primitive  equation  model.  A  similar  forecasting  procedure  is  used  in 
this  model.  The  importance  of  gravity  waves  in  small  scale  convective 
activity  is  discussed  by  Matsumoto  and  Ninomiya  (1969).  Orszag  (1969) 
presented  a  model  for  the  simulation  of  turbulence  which  employs  a 
fast  Fourier  Transform  instead  of  the  Liebmann  Relaxation  Technique. 
Further  refinement  of  this  model  may  be  accomplished  by  the  use  of 
such  a  method. 

Norton  (1970)  discussed  the  development  of  a  three  dimensional  model 
for  application  to  mesoscale  motions  based  on  the  Boussinesq  equations. 
The  model  presented  here  is  a  further  development  of  the  same  model  with 
a  staggered  grid  system  applied  to  minimize  computer  storage  requirements 
The  total  grid  size  is  14  by  14  in  the  horizontal  with  10  levels  in  the 
vertical,  while  the  computational  grid  size  is  11  by  10  by  6.  The  grid 
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interval  in  all  directions  is  1  kilometer.  The  Coriolis  parameter, 
treated  as  a  constant,  is  determined  at  a  mean  latitude  of  35°  North. 
Numerical  computations  are  made  using  the  IBM  OS/360  Computer  of 
the  W.  R.  Church  Computer  Center.  To  conform  to  computational  stability 
a  time  step  of  30  seconds  was  chosen. 


II.  MODEL  DEVELOPMENT  EQUATIONS 


A.  FORECAST  EQUATIONS 

The  pressure  term  in  the  complete  equation  of  motion  can  be  written 
as: 


or 


where 


and 


1 


RT 


-  Vq  p  =  —  v., 
p   3  r   p    3 


p  V3  P  =  0  V3 


6  =  Cr 


Po 


T  < 


e  =  T 


,y< 


Using  these  results  in  the  component  form  of  the  equation  of  motion 
yields: 

&  -  -off  ♦  fv  ♦  vvfu  , 

%  -    .  e  |i  -  fu  +  w|  v  , 


dw        3B. 
dt    '  0  3y 


g  +  vv.,  w  . 


(1) 
(2) 
(3) 


The  perturbation  expressions  can  be  formed  by 
3  =  60(z)  +  6'  (x,y,z,t) 
e  =  e0  +  e(x,y,z,t) 


where 


<< 


and 


<< 


and  the  latter  inequality  corresponds  to  the  Boussinesq  approximation. 
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After  substitution  into  equations  (1),  (2),  and  (3),  the  equations 
of  motion  can  be  written: 

&  ■  "9o  It1  +  fv  +vV3u-  (") 

t  -  -eo  If  - fu  +  -4  v>  v 


with  the  relatively  small   products: 

B  111         a  Ifil         and  ft  ^1 

9  3x       '  9  ay       '  and  e  3Z 

being  neglected.     The  hydrostatic  equation  relates  6     and  e     by: 

36 

_o     =     .  5_ 

9z  e 

o 

and  the  vertical  component  can  now  be  written: 

SI-  -vlf+?  +-3»  <6> 

By  defining  <j>  =  eg'  and  substituting  into  (4),  (5),  and  (6)  and 
expanding  the  total  derivatives,  the  system  of  equations  becomes: 

||  .  _L(U,  .  ||  t  fv  +  vV2  u>  (7) 

3£  -  -L<v>  -  &  - fu  +  -4 »•  <8> 


where 


dt  x   '       3y 

M  ..lm-H  +m  +VV2W,  (9) 

|f    =  -L(e)  ,  (10) 

v3   •   V     =     0  .  (11) 

L(S)     =     -  v3-  (V3S),     S  =  u,v,w,e    . 
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Incompressible  flow  is  assumed,  to  eliminate  sound  waves.  This  system 
has  become  known  as  the  anelastic  equations  and  are  quite  accurate  for 
shallow  convection.  (Ogura  and  Phillips,  1962) 

B.  DETERMINATION  OF  THE  PRESSURE  TERMS 

Equations  (7),  (8),  and  (9)  may  be  combined  in  vector  form  and 
written: 

^=-(V3  •  v)V3  -  v3*  -f(k  x  V3)  +^k  +  v  V^V3     (12) 

By  taking  the  three  dimensional  divergence  and  assuming  the  time  rate  of 
change  to  be  zero,  (12)  becomes: 

fc<'3  •  V  -  0 

=   -V3    C(V3.73)V3]    -    vf* 

-V3.tfOT  x  T3)]  +  |-    ff+VvV^V  (13) 

For  the  horizontal  extent  used  in  this  model,  the  Coriolis  parameter  may 
be  treated  as  a  constant  and  then  (13)  becomes: 

v^  -  -v3-C(Vv3)V3]  -  f0  v3-(k  x  V3)  +  f-  || 


0 

Note  that  when  the  terms  on  the  RHS  are  solved  by  ordinary  differencing 
techniques,  the  value  of  <j>  on  the  LHS  can  be  solved  by  using  the 
Extrapolated  Liebmann  Relaxation  Scheme. 
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III.  FINITE  DIFFERENCING 

Numerical  integration  of  equations  (7)  through  (10)  is  carried  out 
using  finite  difference  methods  in  x,y,z,  and  t.  Extrapolated  Liebmann 
Relaxation  is  used  to  solve  for  the  pressure  term  in  equation  (14).  For 
the  complete  development  of  the  relaxation  scheme,  see  Appendix  A.  The 
technique  applied  to  generate  the  forcing  function  used  in  (14)  is  shown 
in  complete  form  in  Appendix  B. 

To  maintain  consistency  in  finite  differencing  throughout  the  fore- 
cast procedure,  the  differencing  scheme  used  in  equation  (14)  must  be 
equivalent  to  the  type  used  on  equations  (7),  (8),  and  (9).  The  flux 
form  differencing  method  devised  by  Arakawa  (1966)  was  utilized  in  the 
form: 


L(S) 


i,j»k 


(ui+l,j,k  +  ui,j,k 


-(u.  n  .,  +u.  ., 


tvi,j+l,k  +  Vi,j,k 


-(v.  .  ,  ,  +  v.  .  , 
i  ,J-1  ,k    i  ,j,k 


(w.  .  , .,  +  w.  .  . 


(Si+l,J,k  +  Si,j,k 
(Si-l,j,k  +  Si,j,k 


(Si,J+l,k  +  Si,j,k 
(Si,J-l,k  +  Si,j,k 


(Si,J,k+l  +  Si,j,k 


"(wi,j,k-l  +  wi,j,k}  (Si,j,k-l  +  Si,j,k 


/  4ax 


/  4Ay 


/  4az 


This  method  is  used  for  the  u,v,w  and  e  fields,  but  for  the  linear  space 
derivatives  for  <j>,  an  averaging  technique  must  be  employed  due  to  the 
staggered  grid  system.  The  centered  form  of  the  space  derivative  takes 
on  the  form: 
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(—1        =  (S.  .  ,  +  S.  .  .  n  +  S.  .  -,  ,  +  S.  .  ,  ,  J  /  4. 
lax/,,  -,  ,,,    v  i  ,J,k    i,j,k-l    i,j-l,k    i,j-l,k-r  ' 

-(si-l,j,k +  si-l,j,k-l +  si-l,j-l,k +  s1-1,j-1,k-1'/4f4x 

with  similar  forms  for  the  y  and  z  derivatives.  The  primed  subscripts 
refer  to  the  grid  points  one-half  grid  distance  above  or  below  the 
diagonals  connecting  surrounding  velocity  points  in  a  horizontal  plane. 
Centered  time  differencing  ("leapfrog")  is  used  for  all  forecasts  with 
the  exception  of  the  first  time  step  where  a  single  forward  time  step 
is  used. 
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IV.  BOUNDARY  CONDITIONS 

The  boundary  conditions  are  selected  to  enforce  mass  and  energy 
conservation  within  the  confines  of  the  computational  grid.  Figure  1 
illustrates  the  horizontal  grid  domain  and  boundary  placement.  Cyclic 
continuity  of  parameters  is  imposed  by  forcing  periodicity  in  the 
east-west  direction  according  to  the  relations 

S1,J,K  =  SIM-2,J,K 

S2,J,K  =  SIM-1,J,K 

SIM,J,K  =  S3,J,K 

where  S  =  u,v,w,  e.  Since  values  of  S  are  required  at  three  consecutive 
points  in  solving  for  the  pressure  and  divergence  terms,  values  outside 
the  domain  can  be  obtained  through  the  cyclic  conditions.  The  east- 
west  boundary  conditions  on  <f>  are  also  cyclic  and  are  specified  by: 

*1,J,K  =  *IM-2,J,K 

*IM-1,J,K  =  4>2,J,K  . 

The  south  boundary  is  located  halfway  between  velocity  grid  points 
J  =  2  and  J  =  3,  while  the  north  boundary  is  between  the  points 
J  =  JM  -  2  and  J  =  JM  -  1 .  To  prevent  the  flow  of  mass  and  energy 
through  these  boundaries,  the  v-component  is  set  by  the  following 
conditions: 

VI,1,K  =  "VI,4,K 
VI,2,K  =  "VI,3,K 
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1  =  1 
i=JM 


i=IM 


Norlh 


=C 


S  ouih 

-v- 


•  JM-1 


•  JM-2 


•   3 


•  2 


i=  1 
ir  1 


3 


•       •       * 

IM-2     |M-1     i  =  IM 


•  •u.v.w.e 


Figure  1.  Illustration  of  the  grid  system  in  the  horizontal 
direction. 
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VI,JM,K  ~  ~VI,JM-3,K 

VI,JM-1,K  =  "VI,JM-2,K  . 
The  values  of  u,w,  and  e  are  extrapolated  outward  by: 

SI,1,K  =  SI,2,K  =  SI,3,K 

SI,JM,K  =  SI,JM-1,K  =  SI,JM-2,K 

where  S  =  u,w,e.  This  implies  that  there  is  no  advection  through  this 
boundary  and  furthermore  no  heat  exchange  as  discussed  by  Ukaji  and 
Matsuno  (1970). 

The  geostrophic  relationship: 


Wi'.j'.k' 


=  "  fouv,yx 


was  used  to  obtain  the  values  of  <J>  outside  the  north-south  boundaries. 

* 
Where  u.,    .,   .  ■   represents  the  volume  average  of  the  u  velocity  field  at 

1   j  J  »K 

the  boundary.  The  results  reported  here  were  for  the  special  case  of 
f  =  0.  For  this  case,  the  boundary  condition  then  becomes: 


(4. 


=  0  . 
j\k' 


Figure  2  illustrates  the  vertical  grid  domain  and  boundary  placements. 
The  top  and  bottom  boundaries  are  set  in  manner  similar  to  the  north- 
south  boundaries  where  the  points  midway  between  K  =  KM  -  1  and  K  =  KM  -  2 
and  K  =  2  and  K  =  3  are  the  top  and  bottom  boundaries  respectively.  The 
no-flux  condition  is  imposed  on  the  w  field  by: 

WI,J,1   =  -WI,J,4 
WI,J,2  =  "WI,J,3 
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,=  1  -  '=,M 

k=KM 


A 


#  .  •   KM-1 
Top 

\ ■ 

•  •  • KM-2 


#  •  •    KM^3 

•  •  • 


•     3 


Bollom 

— V— 


•     2 


|=  1  2  3  IM-2  IM-1  i=IM 

k=1  ..(J  k=1 


•  -u.v.w.e 


Figure  2.     Illustration  of  the  grid  system  in  the  vertical 
direction. 
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WI,J,KM     ~WI,J,KM-3 


WI,J,KM-1   ~  ~WI,J,KM-2   , 


while  the  e  values  are  set  by  the  relations: 


0I,J,1  "  _6I,J,4 


9I,J,2     "9I,J,3 


9I,J,KM    ~9I,J,KM-3 


9I,J,KM-1  =  '9I,J,KM-2  . 
The  values  of  u  and  v  are  simply  extrapolated  outward  by: 

SI,J,1  =  SI,J,2  =  SI,J,3 

SI,J,KM  =  SI,J,KM-1     SI,J,KM-2 

where  S  =  u,v. 

The  value  of  the  pressure  term  is  obtained  by  the  use  of: 


tel 


The  importance  of  proper  boundary  conditions  cannot  be  over- 
emphasized. This  model  was  found  to  be  extremely  sensitive  to  all  the 
boundary  conditions.  The  complete  development  of  the  boundary  conditions 
is  found  in  Appendix  C. 
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V.   INITIAL  CONDITIONS 

Initialization  of  the  u,v,w  and  e  fields  is  accomplished  by  the  use 
of  linearized  gravity  wave  solutions,  which  are  exact  for  small  amplitude 
gravity  waves  and  are  derived  from  the  following  equations: 

au  _    z± 

at      ax 

3V        3<fr 

at    "  ay 

aw  -    3<t>  .  ge 

at    '  az    e 

o 


and 


v3  •  V  =  0 


M.  =  _  w  — 

3t  3Z 


The  exact  solutions  for  a  stable  wave  propagating  in  the  x  direction  are 

given  by: 

xmW  /  \ 

u     = p y  sin(x  x  -ct)  cos  \xy  cos  mz     , 

X     +  y 

v     = 7p^ — P  cos(x  x  -ct)   sin  \iy  cos  mz     , 

X     +  u 

w     =     W  cos(x  x  -ct)  cos  vy  sin  mz     , 

9     =     -r—    — -  sin(x  x  -ct)  cos  \iy  sin  mz     . 

AC         0  £- 

where  c  represents  the  phase  speed  of  the  stable  gravity  wave  in  m  sec" 


and  is  calculated  from: 


i 


c  =  -?-  \\  -  x2  !  f2  2  (15) 

X  +  y  +  m 
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and  N  is  the  Brunt  Vaisala  frequency 


■Vi 


"  =*'t  H 


The  initial  guess  for  the  $  field  is  chosen  to  be  an  arbitrary  constant 
(2.0)  due  to  the  extremely  efficient  relaxation  technique.  The  set  of 
exact  solutions  with  t  =  0  was  used  to  initialize  the  fields. 
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VI.  RESULTS 

By  introducing  a  stable  gravity  wave  perturbation  into  the  basic 
model,  the  wave  form  can  be  tested  for  proper  phase  speed,  amplitude 
growth  or  decay  and  computational  stability.  The  initial  conditions 
should  result  in  a  wave  that  propagates  in  the  x  direction  alone  with 
no  growth  in  amplitude.  Figures  3a,  3b,  3c,  and  3d  are  representative 
plotted  values  of  the  u,v,w,  and  e  wave  forms  after  time  steps  10, 
20,  30,  and  40,  with  the  initial  value  present  for  comparison.  These 
plots  are  for  the  wave  forms  in  the  x  direction  for  J  =  5  and  K  =  7 
and  show  the  phase  speed  relationship  between  the  various  time  steps. 

The  exact  solution  to  equation  (15)  yields  a  phase  speed  of  5.5  m  sec" 

-1  -3    -1 

for  the  specified  values  W  =  .01  m  sec   and  N  =  4.04  x  10   sec  . 

Phase  speed  calculated  from  the  plotted  values  is  5.4  m  sec"  for  each 
of  the  four  wave  forms.  There  is  negligible  amplitude  growth  or  decay 
in  any  of  the  plots. 

Figures  4a,  4b,  4c,  and  4d  are  the  plotted  values  of  the  wave  forms 
in  the  y  direction  for  I  =  5  and  K  =  7.  The  plots  indicate  there  is 
no  propagation  in  the  y  direction.  Due  to  the  nature  of  the  solutions 
to  the  linearized  equations,  the  wave  forms  should  oscillate  about  the 
zero  value  as  indicated  in  the  stated  figures. 

Plots  of  the  wave  forms  in  the  z  direction  for  I  =  5  and  J  =  5  are 
shown  in  Figures  5a,  5b,  5c,  and  5d  and  are  similar  in  amplitude 
oscillation  to  the  y  direction  with  no  phase  speed  observed. 

The  main  objective  of  this  staggered  grid  model  was  to  reduce  the 
storage  size  and  computer  time  required  for  an  operational  small  scale 
convective  model.  The  previous  full-grid  model  required  502K  bytes  of 
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storage  and  approximately  20  minutes  of  computer  time  to  produce  a  40 
time  step  forecast.  The  forecasts  produced  by  the  staggered  grid 
model  accomplishes  a  40  time  step  prognosis  in  less  than  10  minutes  of 
computer  time  and  requires  only  156K  bytes  of  storage. 

The  results  obtained  after  40  time  steps  compared  to  the  analytical 
solutions  illustrate  the  accuracy  of  the  finite  difference  solutions. 
However,  a  few  time  steps  later,  the  solutions  fail  to  converge.  A 
further  discussion  of  this  difficulty  will  be  given  in  a  later  section. 
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VII.  CONCLUSIONS  AND  SUMMARY 

A  ten-level  staggered  grid  primitive  equation  model  was  applied  to  a 
mesoscale  grid  and  tested  with  an  initial  stable  gravity  wave  perturba- 
tion to  check  the  computational  stability  and  finite  difference  form  of 
the  equations  of  motion.  Comparing  the  results  presented  here  to  a  non- 
staggered  grid  model,  the  computer  time  and  storage  requirements  have 
been  reduced  to  about  1/3  of  the  full  grid  model. 

The  forecasts  produced  by  this  model  compare  quite  closely  with  the 
full  grid  model  but  ceases  to  function  after  47  time  steps.  Up  until 
this  time  step  all  fields  appear  to  be  predicted  properly.  Tests  were 
run  in  double  precision  to  determine  if  truncation  error  was  a  problem 
and  the  same  type  of  results  occured  with  the  model  failing  to  converge 
after  47  time  steps.  In  all  cases  the  points  that  failed  to  converge 
were  near  the  boundaries  which  leads  to  the  conclusion  that  further 
refinements  of  the  boundary  conditions  may  be  necessary.  Experiments 
have  shown  that  a  less  stringent  relaxation  tolerance  produces  less 
exact  $  gradients  which  in  turn  induces  greater  errors  in  the  prognosis 
of  the  other  parameters.  Apparently  there  is  still  some  difficulty  in 
the  program  and  presently  studies  are  being  co. ducted  to  further  refine 
this  model . 

Continued  studies  using  this  model  should  have  a  provision  for 
including  vertical  compressibility.  Additionally,  since  a  primary 
source  of  the  energy  in  a  thunderstorm- type  cloud  is  latent  heat, 
application  to  cumulus  scale  convective  activity  would  be  inadequate 
without  the  inclusion  of  moisture. 
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APPENDIX  A 
EXTRAPOLATED  LIEBMANN  RELAXATION  SCHEME 
FOR  THREE  DIMENSIONAL  APPLICATION 

The  equation  to  be  solved  by  relaxation  is: 

72<f>  -  F  =  0 
where  the  tendency  of  the  divergence  is  zero  and: 


(A-l) 


F  =  -v. 


(V3  ■  v3)V3 


+  foc  + 


3_  M. 

e„  3z 
o 


(A-2) 


with  the  <j>  points  staggered  in  the  vertical  as  well  as  in  the  horizontal, 
(They  are  located  one-half  grid  distance  above  or  below  the  diagonals 
connecting  velocity  points  as  in  Fig.  Al.)  In  other  words,  the  velocity 
and  e  points  are  located  at  the  corners  of  a  cube  and  the  <j>  point  in  the 
center.  For  convenience,  a  new  temporary  indexing  scheme  will  be 
utilized  by  writing: 


i 


=  l 


y   -  J 

k'  =  k 


where  i,j,k  represent  u,v,w  and  e  points  and  i',  j1,  k'  represent  <j> 
points.  Writing  (A-l)  with  this  new  notation  in  finite  difference  form 
at  the  point  i',  j',  k1  and  expanding  yields: 


1 


(4ax) 


'i'+l.j'+l.k'+l   +  2<V+l,j',k'+l   +  *i'+l,j'-l,k'+l   "  2*i',j+l,k,+l 


~4*i\j\k'+l    "  2*i'fj'-l,k'+l  +  ♦I'-l.J'+l.k'+l   +  2*i'-l,j*,k'+1 
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Figure  Al .  Three  dimensional  schematic  of  the  location  of 
Various  fields. 
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(4Ay)' 


+V-l,j'-l,k'+l  +  2*i  '+1  ,j'+l  ,k'  +  4i,+lj,,k'  +  2*i  '+1  .j  '-1  »k' 

"4*i',j+l,k'  "8*i',j',k'  "  4*i',jM,k'  +  2*i'-l,j'+1,k' 
44*1,-l,j,.k'  +  ^I'-l.J'-l.k'  +  ♦l'+lfJ,+l.k,-l  +  ^i'+l.j'.k'-l 
^i'+l.j'-l.k'-l  "  ^i'.j'+l.k'-l  "  4*i',j,ik,-l  '  2*1" .J'-l.k'-l 

+4,i'-l,j'+l,k'-l  +  2*i'-l,j',k'-l  +  ♦l'-l,j,-lik,-l  J  + 

♦I'+l.j'+l.k'+l  "  2*i'+l,j,k'+l  +  *i'+l,j'-l,k'+l  +  2<V,j'+l,k'+l 
"^I'.j'.k'+l  +  2<V,j'-l,k'+l  +  ♦l,-ltJ,+ltk,+l  "  ^I'-l.J'.k'+l 
+*  I'-l.j'-l.k'+l  +  ^I'+i.j'+l.k'  "  4*i'+i,j',k'  +  2<J,i'+l,j'-l,k' 
+4*i',j'+l,k'  "  8*i',j\k'  +  4*i',j'-l,k'  +  2*i'-l,j'+l,k' 
"^I'-l.j'.k'  +  2*i,-1,J,-l,k'  +  ♦ll+l,J'+l,k,-l  "  ^I'+l.j'.k'-l 
+V+l,j'-l,k'-l  +  ^I'.j'+l.k'-l  ■  ^I'.J'.k'-l  +  2*i\j'-l,k'-l 


+($>i  i    i     -;  ixi    m    i    ~    2d)_. ■* 


(4az)' 


T-l.j'+l.k'-l       ^i'-l.j'.k'-l       T-l.j'-l.k'-l 

♦I'+l.j'+l.k'+l  +  ^i'+l.j'k'+l  +  *i,+lsj,-l,k,+l  +  2*i '  ,j  '+1  ,k'+l 
+4*i' .J'.k'+l  +  2<V,j'-l,k'+l  +  ♦I'-l.j'+l.k'+l  +  2*i'-l,j',k'+l 
^V-l.j'-l.k'+l  ~  2V+l,j'+l,k'  "  4*i'+l,j',k'  "  ^i'+l.j'-l.k' 


-4*i',j'+l,k'  "  8*i\j',k'   "  4*i',j'-19k'  "  2*i'-l,j'+l 


k' 


'i'-l.j'.k"  "  2*i'-i,j'-i,k'  +  *1,+l,j,+l,k'-l  +  ^V+l.j'.k'-l 
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+*r+i,j'-i,k'-i  +  2*i',j'+i,k'-i  +  4*i',i',k'-i  +  2*i',j'-i,k'-i 


^i'-l.j'+l.k'-l  +  ^i'-l.j'.k'-l  +  V-l.j'-l.k'-l 


-Fi',j',k'  =  ° 


(A-3) 


Rearranging  terms  and  clearing  fractions  leads  to 


(4Ay)2  (4az)2 


+(4ax)2  (4az)2 


+(4ax)2  (4Ay)2 


B  -  84,.ljjIjk, 


C  -  8*ilJljkl 


D-8,iljjl)k, 


-(4ax)2  (4Ay)2  (4az)2  F.,^,^,  =  0 


(A-4) 


where  B,  C  and  D  represent  the  finite  difference  forms  in  equation  (A-3) 
with  the  term  84? . ,  .,  ■  ,  subtracted  off.  With  further  rearranging  of 
terms  ^-4)  becomes: 


♦v.j'.k'  =  A 


(4Ay)2  (4az)2  B  +  (4ax)2  (4az)2  C 


+  (4ax)2   (4Ay)2  D  -    (4ax)2   (4Ay)2   (4az)2  F . ,    .,    . 

1    » J    » i> 


where: 


A  = 


8  (4Ay)2  (4az)2  +  (4ax)2   (4az)2  +  (4ax)2   (4Ay)2 
old 


Adding  and  subtracting  <|> . ,    .,    .  ,  yields 


* 


new 


?ld..,   fc,«A 


i'.j'.k'       M\j\k 


(4Ay)2  (4az)2  B  +  (4ax)2  (4az)2  C  + 


+  (4ax)2   (4Ay)2  D  -    (4ax)2   (4Ay)2   (4az)2  F . ,    • ,    ,, 

I         5  J        >  *• 


(A- 5) 
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Where  the  over-relaxation  coefficient  a  is  introduced.   (A-5)then 
becomes: 


new 


old 


*i\j',k'  =  H  -  a)  V,j',k'  +a  A 


(4Ay)2  (4az)2  B 


+  (4ax)2  (4az)2  C  +  (4ax)2  (4Ay)2  D 
-  (4ax)2  (4Ay)2  (4az)2  F.,^.,^,1 


(A-6) 


This  was  the  form  used  for  relaxation  in  the  model.  Many  tests  were  run 
to  determine  the  most  efficient  relaxation  coefficient  for  this  type  of 
finite  difference  scheme.  The  value  of  a  which  resulted  in  the  most 
efficient  relaxation  was  1.05.  Note  that  if  ax  =  Ay  =  az,  equation  (A-6) 
reduces  to  the  more  familiar  form: 


new 


*i',j\k'  -  (1  *  a)  ♦I'.J'.k'  +  24 


B  +  C  +  D 


-  (4azT  F 


i'.j'.k' 


The  stencil  of  the  method  of  the  *  operator  in  the  three  dimension- 
al staggered  grid,  assuming  that  the  three  grid  lengths  are  equal,  is 
shown  in  Fig.  A2. 
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K-1 


Figure  A2.  f2  Operator  in  three  dimensional  form  assuming 


AX  =  Ay  =  AZ, 
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APPENDIX  B 
DEVELOPMENT  OF  THE  FORCING  FUNCTION 


The  formation  of  the  forcing  function  used  in  the  relaxation  scheme 
is  a  long  process.  Consistency  between  the  finite  difference  methods 
used  in  the  prognostic  equations  and  in  the  relaxation  technique  must 
be  maintained.  The  computation  of  the  forcing  function  given  in 
equation  (A-2)  in  finite  difference  form  follows.  In  order  to  determine 
this  forcing  function  in  a  consistent  manner,  the  forcing  function  must 
first  be  calculated  at  an  i,j,k  point  and  then  volume  averaged  to 
have  the  function  at  an  i',j',k'  point.  Each  of  the  three  directional 
components  are  computed  separately  and  then  summed  to  put  the  complete 
forcing  function  at  the  desired  point.  For  example,  in  order  to  form 
Fgi  ri  n  (the  forcing  function  at  the  point  i '  =  5,  j '  =  5,  k'  =  5) 
the  following  steps  are  required: 


FX 

J  1 

~  (  4 

1 

"  4 

Fy 
r5',5',5' 

(l 
(4 

1 

"   4 

Fz 
r5\5',5' 

(l 
"  (4   . 

1 
"  4 

G(U)6.6.6  +  G(U)6,5,6  +  G(U)6.6.5  +  GH.5.5 


G(UW  +  G(u)5,5,6  +  G(u)5,6,5  +  G<U>5,5,5 
G<VW  +  GW5.6,6  +  G^6,6,5  +  ^5,6.5 

G<V>6,5,6+G<V>5,5,6+G<V>6,5,5+G<V>5,5,5 


BM6.6.6  *  G<w>6.5.6  \GM5.6.6  *  GM5.5.6 


G<w)6,6,5+G(w>6,5,5+G(w>5,6,5+G<w>5,5,5 


I  i- 

(     AX 


1_ 

(Ay 


h 

|AZ 


F  ,   ,   ,   =  Fx,   ,   ,   +  Fy,   ,   ,   +  FZ, 

•J      5  O      j  u  J   5  U   5  0  0)0)0  0)0)0 
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where 

G(u)  =  -  L(u)  +  fv  +  vv^  u  , 

2 
G(v)  =  -  L(v)  -  fu  +  vv3  v  , 

6(w)  =  -  L(w)  +  f-^  +  vv*  w  . 

eo    3 

and  the  G  functions  are  calculated  on  the  i,j,k  grid 


44 


APPENDIX  C 
COMPLETE  DEVELOPMENT  OF  THE  BOUNDARY  CONDITIONS 

For  the  boundary  described  by  the  plane  of  j '  =  2*  the  y  component 
of  the  equation  of  motion  may  be  written 

H       =0 

3ti',2,k' 


which  implies 


3V  3V 

"l,3.k     "  3ti,2,k  {C-]) 


Further  expansion  of  each  of  the  terms  leads  to: 


—     =  -  LM      -  2£  _  f  u  (C-2 

9ti,3,k       ;i»3,k  WU3ik       Vi,3,k  lt  c 


and 


-  n  =   L(v)i  2  k  +  f£     +  Vi  2  k  (c"3) 
9ti,2,k      1'^»K   3yi ,2,k    °  1'*»K 

By  introducing  the  no-flux  conditions  on  v  at  the  lateral  walls  from 
(C-l),  (C-2)  and  (C-3)  yield: 

-  $&■  -fu      =  $±  +  f  u  (C-41 
9yi,3,k    °i'3'k   9yi,2,k    °  i'2'k 

Expansion  of  (C-4)   into  finite  differences  and  averaging  the  u 

field  to  be  consistent  with  the  4>  field  yields: 

f 

♦I'.l.k'    =  *V,3#   +  HP     (ui,2,k  +  Vl,2,k  +  ui,2,k+l 

♦"WW1 
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with  the  no-flux  conditions  being: 

vi,2,k     =     ~vi,3,k 

Vi']'k  =  "vi,4,k 
A  similar  approach  is  taken  at  the  northern  boundary  with  the  following 
conditions: 

vi,Jl,k  =  ~vi,J2,k 
vi,JM,k  =  "vi,J35k 

*i',Jl,k*  =  *i'  ,J3,k'  "   2  (ui,J2,k  +  ui+l,J2,k 

+ui,J2,k+l  +  ui+l,J2,k+l  *   * 

The  case  considered  here  was  for  f  =  0. 

For  the  boundary  described  by  k1  =  2',  the  z  component  of  the 
equation  of  motion  may  be  written: 

3V,2',k' 

which  implies: 

3W  9W 

3ti,3,k      3ti ,2,k 
Further  expansion  of  each  of  the  terms  in  (C-5)  yields: 


(C-5) 


il     =  "  L(w)i  3  k  "  if     +  f-  e-  3  k  (C-6) 

9ti,3,k       1,J'K    3zi,3,k   9o   1,J,K 

il     =  +  L(w)i  2  k  +  if     "  f-    ei  2  k  •  (C'7) 

9Ii,2,k       ud*K       3zi,2,k   eo   l5^'k 
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By  introducing  the  no-flux  conditions  on  w  at  the  top  and  bottom 
boundaries  from  (C-5),  (C-6)  and  (C-7)  they  yield: 

3zi,3,k   9o  i'3'k   3zi,2,k   6o   1»2'k  {       } 

By  setting  e  at  the  boundaries  as  previously  discussed,  (C-8)  is 
satisfied  if: 

where  the  9  boundary  conditions  are: 

8  •  •  /->   =   -  9  .  •  o 

i*J.2       i,j,3 

6i,j,l  =  "  9i,j,4  . 

A  similar  approach  is  taken  at  the  top  boundary  with  the  following 
conditions: 

9i\J,Kl     =     "9i,j,K2 


and 


9i,J,KM     =     ~9i,j,K3 


V.j'.KT      =     V.JMQ' 
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